Transient behavior of full counting statistics in thermal transport 
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The generating function of energy counting statistics is derived for phononic junction systems. 
It is expressed in terms of the contour-ordered self-energy of the lead with shifted arguments, 
E a (t,t') = T*l(t + hx(r),r' + Tix(t')) — £i,(t,t'), where Si(t,t') is the usual contour-ordered 
self-energy of the left lead. The cumulants of the energy transferred in a given time tu from the 
lead to the center is obtained by taking derivatives. A transient result of the first four cumulants of 
a graphene junction is presented. It is found that measurements cause the energy to flow into the 
lead. 
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Phonon transport in the ballistic quantum regime pos- 
sesses special features, such as the quantized universal 
thermal conductance 0, 0] and wave-like coherent trans- 
port described by a Landauer-like formula @, 0] . A typ- 
ical set-up of such a system consists of two infinite heat 
baths maintained at different temperatures with a finite 
junction part forming the scattering region. The focus 
in the last decade has been on steady-state thermal cur- 
rents. Since the heat baths are stochastic in nature, it is 
natural to ask a statistical question: what is the distribu- 
tion of the energy Q transferred in a given time tu- Such 
questions have been raised in electron transport, where 
it is known as the full counting statistics. Levitov and 
Lesovik presented their celebrated formula which forms 
the definite answer to the question [5j. Many works fol- 
lowed in electronic transport [1,0]- The electron counting 
statistics has been experimentally measured in quantum 
dot systems [8[ . No such measurements have been carried 
out for thermal transport, but it is potentially possible, 
e.g., in a nano-resonator system. 

Saito and Dhar [j| treated the full counting statistics 
for heat transport in a ID chain. Such inquiries also 
have deep connections with the nonequilibrium fluctu- 
ation theorems [l(J. The result obtained by Saito and 
Dhar was only for the long-time limit. In this paper, we 
present a formulation based on two-time measurements, 
treating the transient behavior and long-time limit on an 
equal footing. A central result of our derivation is that 
the generating function can be concisely expressed by 
the contour-ordered self-energies of the lead, making con- 
tact with the nonequilibrium Green's function (NEGF) 
method Q of quantum transport. A more general expres- 
sion for the long-time limit of a general junction system 
with any number of degrees of freedom is also derived, 
and numerical results for the transient behavior of the 
first few cumulants of a graphene junction are presented. 

We consider initially decoupled harmonic systems de- 
scribed by the Hamiltonians 



for the left and right leads and a central region. The 
leads are assumed semi-infinite while the center has a 
finite number of degrees of freedom. Masses are absorbed 
by defining u — ^pmx. u a and p a are column vectors of 
coordinates and momenta. K a is the spring constant 
matrix of region a. Couplings of the center region with 
the leads are turned on either adiabatically from time t = 
— oo j or switched on abruptly at t = 0. The interaction 
term takes the form H int = u£V LC uc + UgV RC uc- The 
total Hamiltonian is H = Hr, + He + Hr + H lnt . 

Focusing on the left lead, we define the energy current 
operator by the rate of decrease of energy of the lead (in 
the Heisenberg picture) as 



dH L {t) = %_ 



= = -[H L ,H H ] = PL (t) T V LC u c (t), 



dt 

(2) 

where Hjj is the Hamiltonian in the Heisenberg picture. 
We define the 'heat' operator as 

Q = [ I(t') dt' = H L - U{0, t)H L U{t, 0), (3) 
Jo 

where Hl [= Hl(0)] is the Schrodinger operator of the 
free left lead, and U(t, t') is the evolution operator under 
H(t). U satisfies the Schrodinger equation 



ih^p-=H(t)U(t,t'). 



(4) 



What we would like to calculate is the moments of the 
heat energy transferred in a given time t. To this end, we 
look at the generating function of the moments instead. 
Since Q is a quantum operator, there are subtleties as to 
how exactly this generating function should be denned. 
Naively, we may use (e 1 ^). But this definition fails the 
fundamental requirement of positive definiteness of the 
probability distribution, 



P{Q) 



2tt' 



(5) 



Ha = -^P a Pa 



1 



u T a K a u a , a = L,R,C, (1) 



for a classical quantity Q. The correct definition is 

Z= ( e ^ e -^ M)', (6) 



2 



based on measurements at time and t where each time 
a measurement of the energy of the left lead is carried 
out, the wavefunction collapses into the eigenstate of the 
operator H L . Thus, to take care of this process, the 
average is defined by 



>'=Tr[£P o p(0)P o - 



(7) 



where P a is the projector onto the eigenstate of Hl with 
eigenvalue a. p(0) is the steady-state density operator 
obtained by adiabatically evolving from a product state 
at t = — oo to t = 0. 

To calculate the generating function Z, we use the 
following strategies. First, the projector is repre- 
sented by Fourier transform, P a = 5 (a — Hl) — 



I: 



-i\(a-H L ) 



d\/{2-K). Then, the products of the ex- 



ponential factors in Z, combined with the exponential 
factors in the projectors, are written in terms of an evo- 
lution operator U x (t,t') of an effective Hamiltonian with 
a parameter x, given 



z(0 



, e iZH L /2 e -iiH L (t) e tiH L /2y 
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The proportionality constant will be fixed later by the 
condition Z(0) = 1. The evolution operator U x is associ- 
ated with the Hamiltonian 



H x (t) = e ixHL H(t)t 



— ixHr 



H(t) + (u L (hx)-u L ) V LC u c , (9) 



where Uhi'hx) — e txHL ULe~ lxHL is the free left lead 
"Heisenberg" evolution to time t = fix. We can give 
a more explicit form for the Hamiltonian, 



H x (t) = H{t) + [u T L C{x) + plS(x)] uc, (10) 



where 



tLC 



C(x) = (cos{Tix^/Kl) - \)V L 
S{x) = (1/ \J~Kl) sin(hxy/K^)V LC . 



(11) 

(12) 



Next, we represent U x using path integrals. The la- 
grangians associated with the path integrals are (ignoring 
the right lead for the moment): 



r 1-2 1 T rsL 

C c = \ul-\u T c (KC~S T S)u Cl 
Clc = -ulSu c - u T L (V LC + C)u c 



(13) 

(14) 
(15) 



Following Feynman and Vernon 12 1, we can eliminate 
the leads by performing gaussian integrals. Since the cou- 
pling to the center is linear, the result will be a quadratic 



form in the exponential, i.e., another gaussian. The in- 
fluence functional is given by 

I[u c (t)\ = [ V[u L ]p L (-oo)e* J ' d ^+^c) 



Tr 



Vi(t) 



Z L 

K ff drdr'ul(T)n(ry)uc(T') 
1 



(16) 



%x(t))V 



LC 



-UcS+Suc (17) 



In the above expressions, the contour function uc( T ) is 
not a dynamical variable but only a parametric function. 
T c is the contour order operator. Note that Vi is the 
interaction picture operator with respect to Hl , as a re- 
sult, e ltHL / h u L {fix)e~" ltHL / Tl = u L {t + hx). We define the 
contour function x(t) as whenever t < or t > tu- 
Otherwise it is x + (t) — — £/2 — A on the upper branch, 
and x~ (t) = £/2 — A on the lower branch. The important 
influence functional self-energy on the contour is 



H(t,t') 



yA 
^L 



S t SS(t, r'), 



(18) 



LC 



S A + S L = V CL g L {r + Tix(t),t' + Tix(t'))V 

= Z L (T + hx(T),T' + Kx(t')), (19) 

where £_l is the usual lead contour self-energy, S is the 
Dirac delta function defined on the contour. Equation 
([T9"l) is the most important equation defining the self- 
energy of the problem. The generating function Z can 
be expressed in terms of the usual Green's function G = 
G cc of the central region and this particular self-energy. 
The self-energy T, A is obtained from the lead self-energy 
by appropriately shifting the contour time arguments 
and taking a difference. With this result, infinite degrees 
of freedom (due to the semi-infinite nature of the leads) 
reduce to finite degrees of freedom. 

The generating function is obtained by another gaus- 
sian integral, given 

Z((,X) = I V[u c }p c (-^)e {i/h) I dTCc I[u c } 



= J 2?[u c ]p c (-oo)e« Sc - 
ocdetp)- 1 / 2 , 



(20) 



where 



S cS = l -JdrJ dr'u c (T)D(T,r')u c (T'), (21) 



D(t,t') = --^5(t,t')-K c 6(t,t') 



<9t 2 ' 

-£(t,t')-£ A (t,t') 



(22) 



where £ = + We define the Green's function G 
by DqG = 1, or more precisely 



D q {t,t")G{t" ,r')dr" = 5(t,t>). (23) 
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In the above formula for Z, we imagine that the differ- 
ential operator (integral operator) D and Dq 1 are rep- 
resented as matrices indexed by space j and contour 
time r. We can make a systematic expansion in term 
of Ti A by noting the following formulae for matrices, 

det(M) = e TrlnM , and ln(l - y) = -££=i Using 
this, we can write 



lnZ(C) 



oo 



hm y — Tr 

k=l 



U,r) 



(GE 



(24) 



This formula is the central result of this paper. The 
expression is valid for any transient time tM embedded 
in the self-energy Yi A . The notation Tru^ means trace 
both in space j and contour time r, i.e., integrating over 
the Keldysh contour. The projection to the eigenstates 
of Hl results in an integration over A. Since the range 
of the integration is from — oo to +oo, and the two- 
parameter generating function Z(£, A) approaches a con- 
stant as |A| — > oo, the value of the integral is dominated 
by the value at infinity. Our choice of the proportionality 
factor satisfies the required condition of Z(0) = 1. 

For NEGF notations and relations among Green's 
functions, we refer to Ref. Q. It is more convenient 
to work with a Keldysh rotation for the contour or- 
dered functions, keeping Tr(^4i? • • • C) invariant. For any 
A aa (t,t'), with a, a' — ± for branch indices, the effect 
of the Keldysh rotation is to change to 



A = 



A r A K 

A K A a 



A t -A t -A < + A>, A 1 
A t + A*-A < -A>, A* 



(25) 



A* 
A* 



A< + A> 
A< - A> 



We should view the above as defining the quantities A r , 
A a , A K , and A K . For the usual Green's function G we 
get 



G 



G r G K 
G a 



(26) 



The G K component is due to the standard relation 
among the Green's functions. But the K component is 
nonzero for E" 4 . 

In the long-time limit, translational invariance is re- 
stored for the self-energies. Convolution in time domain 
simply becomes multiplication in the frequency domain. 
The shifts given to the arguments in become indepen- 
dent of time t, only depend on the branches. We have 



y t _ yt 



A = 0, 



E2(t) 



(27) 
(28) 
(29) 



Fourier transforming the lesser and greater self-energies, 



we obtain 



we 



Me 



1 



M 



f) . We can now compute the matrix 




FIG. 1. The structure of a graphene junction with 6 degrees 
of freedom with two carbon atoms as the center. 



product GE A . Finally, the generating function for large 
t M is 

/ + OC 7 
pirln(l-GE A ) 

= -*m J ^^Indc^l - G r T L G a T R [(e^-l)f L 

+ (e-^-l)f R + (e^ + e- l ^-2)f L f R ]y (30) 

where G, E A , and T a = i(E„ — E°) are in the frequency 
domain and f a = l/( e ^ Ruj -l), /3 a = l/(k B T a ), is the 
Bose distribution function. This result generalizes that 
of Saito and Dhar Q . It satisfies the steady-state fluctu- 
ation theorem 0, Z(£) = Z(-£ + i(/3 R - /3 L )). 

The long-time result does not depend on how the ini- 
tial states are prepared before measurement. This is 
not the case for transience. The generating function, 
Eq. (|24|) . is for the case where the system is prepared 
in a steady state. A measurement at time disturbs the 
system, and similarly at time tM- Instead of a steady 
state, we can also prepare the system in a product state, 
p(—oo) oc exp(— /3 a H a ). This means that the cou- 
pling iJ; n t is switched on suddenly. Then the projector 
P a commutes with the density matrix with no effect on 
p{— oo). This simplifies the problem. We use the Feyn- 
man diagrammatic technique to obtain the result. Omit- 
ting the details, we have 



InZo 



rTr 



ln(I-GnE y 



(31) 



This expression looks formally the same as before except 
that Go satisfies a Dyson equation defined on the contour 
from to tM and back, while G is defined on the Keldysh 
contour from — oo to tM- 

G (t,t')= 9c (t,t') (32) 

dn dr% gc (r, n ) E (n , r 2 ) G (r 2 , r' ) , 



where gc is the contour ordered Green's function of the 
isolated center. 

We now present some numerical results. Fig. I is the 
structure of our graphene junction system. The center 
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FIG. 2. The cumulants ((Q n >> for n = 1, 2, 3 and 4. The 
curves are for the product initial state; the circles are for 
steady-state initial state. The dotted line is for the classical 
limit (h — > keeping A finite) for the steady-state initial con- 
dition. The temperature of the left lead is 330 K and that 
of the right lead is 270 K. For the product initial state, the 
center temperature is 300 K. 



steady state, the effect of measurement is always to feed 
energy into the measured (left) lead, even if the tem- 
perature of the left lead is lower than that of the right 
lead. If the system were classical, the measurement can- 
not disturb the system. We should expect the current 
to be constant once the steady state is established. The 
nonlinear tu dependence observed here in (Q) is funda- 
mentally quantum-mechanical in origin. 

In summary, the generating function for phononic junc- 
tion systems is obtained, which can be written compactly 
using Green's function as InZ = -(l/2)Trln(l - GS A ). 
A central quantity is the self-energy Y, A which is ex- 
pressed in terms of the usual lead self-energy with shifted 
arguments. This is a very general result valid for steady- 
state initial states or product initial states in a two-time 
measurement. Numerical results for a graphene junction 
system are presented. An intriguing feature is that a 
measurement, even in the steady state, causes energy to 
flow into the leads. We hope that such robust features 
can be verified experimentally. 

This work is supported in part by a URC research grant 
R-144-000-257-112 of National University of Singapore. 



region consists of two atoms, while the two leads are 
symmetrically arranged as strips (with periodic bound- 
ary conditions in the vertical direction). We obtained 
the force constants using the second generation Brenner 
potential. To compute the transient results, we need to 
perform convolution integrations in the time or frequency 
domain many times. It is handled by treating the con- 
volutions as matrix multiplications. Then the expression 
of the derivatives, ((Q n )) — d n \nZ/d(i£ ) ) n , is calculated. 
Note that the £ dependence only enters through T, A . We 
also note that a power series in GT, A terminates after n 
terms for ((Q n )) for the product-state initial condition, 
but it is an infinite series for the steady-state case. The 
computational effort required for convergence is huge for 
the graphene junction. We also obtained the result for 
ID chain which will be presented elsewhere. 

Fig. 2 shows the first four cumulants. The first cumu- 
lant, which is also the first moment, is the total amount 
of energy entering the center from the left lead during 
time to tjv/- Its derivative gives the current. Such tran- 
sient currents have been calculated [H| for the product 
initial states for ID chains. The second cumulant gives 
the variance of Q. The higher order cumulants are small 
but not zero, thus the distribution of Q is not gaussian. 
For large times, all the cumulants become linear in tMi 
and are in agreement to the long-time prediction. 

One striking feature of the results is that the product 
initial state and the steady-state initial state results be- 
have qualitatively the same. The heat transferred, (Q), 
starts from and goes down to negative values. This 
means whether we start from a decoupled system or a 
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